Real-time Monte-Carlo simulations for dissipative tight-binding systems 

and time local master equations 
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The numerically exact path integral Monte Carlo approach for the real-time evolution of dissi- 
pative quantum systems (PIMC), particularly suited for systems with discrete configuration space 
(tight-binding systems), is extended to treat spatially continuous and correlated many-body sys- 
tems. This way, one has to consider generalized tight-binding lattices with either non-equidistant 
spacing or in higher dimensions, which in turn allows to analyze to what extent Markovian master 
equations can be applied beyond the usually studied spin-boson type of models. 
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I. INTRODUCTION 

Quantum Brownian motion is much more involved 
than its classical analog since in general tractable equa- 
tions of motion do not exist pj. Progress can be made 
in two limiting ranges, namely, in the realm of weak dis- 
sipation and in the opposite one of strong friction. In 
the former case a perturbative treatment has led to a va- 
riety of Markovian weak-coupling master equations 0, 
among them the famous Lindblad yj and the Redfield |4| 
equations. Successful applications include quantum op- 
tical systems, decoherence for solid-state based quantum 
bits and nonadiabatic dynamics in molecular systems, to 
name but a few. Strong friction has been explored only 
recently [fj with a growing amount of research since then 
[HQ- There, the quantum Smoluchowski equation, a sort 
of Markovian master equation as well, allows to investi- 
gate condensed phase dynamics at lower temperatures 
e.g. in soft matter and mesoscopic systems |8j. 

A formally exact description of open quantum sys- 
tems valid for all temperatures and damping strengths 
is provided by the path integral approach initiated by 
the work of Feynman and Vernon [9j and developed in 
detail in the 1980s Q, 0, 0]. The approach has been 
utilized in numerous applications especially in condensed 
phase systems, e.g. to reveal the non-exponential decay 
of low temperature correlation functions, and has fur- 
ther been exploited to consistently derive the Markovian 
master equations mentioned above 0, 0] . A real chal- 
lenge, however, has been to evaluate the formally exact 
expression for the reduced density matrix in parameter 
regions where analytical progress and perturbative sim- 
plifications are prohibitive. With the increasing com- 
plexity of designed and controllable quantum systems on 
the nanoscale, e.g. for quantum information processing or 
molecular electronics, the issue of efficient numerical pro- 
cedures in the real-time domain becomes a very crucial 
one. Indeed, important achievements have been made in 
the last decade with the development of advanced path 
integral Monte-Carlo methods (PIMC) the quasi- 

adiabatic propagator scheme (QUAPI) |l4j| . stochastic 
Schrodinger equations 0, and basis set-methods (l6|. 

In this context a certain class of systems has been ex- 



tensively studied, namely, systems with discrete configu- 
ration space, also coined tight-binding systems (TBSs). 
For these systems quantum diffusion takes place on a lat- 
tice, where the sites are coupled by tunneling amplitudes, 
an important case being the restriction to nearest neigh- 
bor coupling. The simplest example is the well-known 
spin-boson model p"?| with applications from condensed 
matter physics to electron transfer reactions. The fun- 
damental role of TBSs follows from the diversity of re- 
alizations in physics and chemistry 0. Transport prop- 
erties in general multistable systems at sufficiently low 
temperatures can be described within TBS models lead- 
ing to relations to the Kondo problem and the Luttinger 
liquid model. Further, TBSs serve as archetypical mod- 
els to study quantum phase transitions in correlated sys- 
tems, as e.g. for Hubbard type of models. Remarkably, 
even a large class of continuous systems can be mapped 
exactly onto TBSs by means of duality transformations 
[l8j . Based on the path integral representation exact 
non-Markovian master equations for TBSs [lj have been 
derived, which provide the starting point for perturba- 
tive treatments such as the non-interacting blip approx- 
imation (NIBA) and reduce to Markovian ones even for 
moderate dissipation and low temperatures |19|. i.e. far 
from the limiting ranges addressed above. 

The PIMC approach is particularly suited to treat the 
dissipative real-time evolution of TBSs numerically ex- 
actly. While former applications were restricted mainly 
to two and three state models [2(j, recently, we substan- 
tially im prov ed the approach to apply to more complex 
systems |2ll |22| . such as single charge transfer across 
long one-dimensional molecular chains including impuri- 
ties and external driving fields. Moreover, the simulation 
time range could be extended to cover basically all rel- 
evant time scales of the dissipative dynamics. We have 
been thus in a position to directly access the range of 
validity of Markovian master equations for TBSs, which 
in turn are extremely helpful to reveal the relevant phys- 
ical processes behind the numerical data. The goal of 
this paper is now twofold: On the one hand we push 
the PIMC procedure even further and present first re- 
sults for the dissipative dynamics of spatially continu- 
ous and of correlated many-body systems; on the other 
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hand, we give arguments to what extent Markovian mas- 
ter equations can be used in these more involved situa- 
tions. We will see that this latter question will directly 
lead us to consider generalized tight-binding lattices with 
non-equidistant spacing between sites or in higher dimen- 
sions. 

The article is organized as follows. We start in Sec. [D] 
with a brief summary of the path integral representation 
for open quantum systems and continue in Sec. lIIII to col- 
lect the main ideas of the PIMC scheme. Sec. I1VI deals 
with known results for non-Markovian and Markovian 
master equations for TBSs, the applicability of which for 
one-dimensional chains is illustrated. Then, in Sec. 
spatially continuous systems are discussed, before we 
come to the correlated many-body dynamics in Sec. IVII 
At the end some conclusions are given. 



1 r , , cosh[w(ft/3/2 - it)] 
dujJ{uj) sinh(ft/? w /2) 



(•5) 



where (3 = l/k B T. 

In the sequel we focus on systems evolving in a dis- 
cretized configuration space with respect to the pointer 
variable q and thus consider the population P{qf, t) of a 
"lattice site" qj determined by the diagonal part of the 
reduced density matrix, i.e., 



P(q f ,t)=Tr{\q f )(q f \p(t)} , 



(6) 



which is normalized J dqP(q,t) — 1. Accordingly, the 
initial density matrix of the total compound W(0) in (J3J 
is taken to be 



W(0) 



Zr% 



[Qi\ e 



-/3(H R - qi fi£) 



(7) 



II. DYNAMICS OF DISSIPATIVE QUANTUM 
SYSTEMS 

The standard approach for the inclusion of dissipa- 
tion into a quantum mechanical formulation starts from 
a system+reservoir model 



H — Ha + Hrt + Hi 



(1) 



with a system part H§, an environmental part Hr, and 
a system-bath interaction Hi. The reservoir (heat bath) 
is mimicked by a quasi-continuum of harmonic oscillators 
bilinearly coupled to the system: 



h r + h i = J2 



pi 
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1 2 



X a 



c a q 



(2) 

where q denotes a system operator corresponding to a 
one-dimensional degree of freedom. Dissipation appears 
when one considers the reduced dynamics by properly 
eliminating the bath degrees of freedom, i.e., 

p{t) = Tr fl {e~ lHt/n W(0) e lHt/h } (3) 

with an initial density matrix W(0) of the total system. 
It turns out that for the reduced dynamics the environ- 
mental parameters enter only via the spectral density 

J ( w ) = irEA^-^)' ( 4 ) 



m a uj a 



which effectively becomes a continuous function of uj for 
a condensed-phase environment. Note that in the above 
definition of the spectral density we have included a fac- 
tor a 2 /h with a being a proper length scale which is 
convenient for the treatment of TBSs. The Gaussian 
statistics of the isolated environment is determined by 
the complex- valued bath autocorrelation function which 
for real time t reads 



with the partition function of the isolated reservoir Zr. 
The bath is equilibrated according to a localized initial 
state of the system on a lattice site qi, where [i£ — 
57J Q c a X a so that e.g. for electron transfer in a polar en- 
vironment, (i is the electronic dipol moment and £ the 
collective dipol moment of the bath. Generalizations to 
delocalized initial states for the system are straightfor- 
ward [13. 

The path integral representation provides a formally 
exact expression for the reduced dynamics and is thus 
the starting point for a numerically exact Monte-Carlo 
(MC) scheme. Along the lines sketched above, the bath 
degrees of freedom are eliminated exactly to arrive at the 
reduced dynamics. As shown in Ref. 0, one thus obtains 
for Eq. © 



jvs S s( t )iSf exp |^5s[s] - $[s] 



(8) 



Here the path integration runs over closed paths s(t) con- 
necting s(0) = Si with s(t) = Sf along the real-time con- 
tour t E — > t — > 0, which combines the forward and 
backward paths s(t') and s'(t'), respectively. Further- 
more, Ss[s] denotes the action of the free system. The 
influence of the traced-out bath is completely encoded in 
the Feynman- Vernon influence functional $[s] 9] 

*[g, q'] = fdt' f dt"W) - q'(t')} [L(t' - t")q{t") 



o Jo 



where 



- L*(t' -t")q'(t")] 
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(9) 



(10) 



The influence functional introduces long-ranged nonlocal 
interactions among the system paths so that in general 
an explicit evaluation of the remaining path integral in 
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Eq. (jSJ is possible only numerically. In this situation 
the PIMC method has been proven as a very promis- 
ing approach to obtain numerically exact results even 
in regions of parameter space where other approximate 
methods fail. 



III. PIMC SIMULATION METHOD 

A prerequisite for an efficient numerical algorithm is 
an appropriate discretization of time and configuration 
space. The latter one is intrinsically given for multi- 
stable systems in the tight binding limit, where the rel- 
evant states are strongly localized in position, only very 
weakly coupled by tunneling, and energetically well sep- 
arated from the rest of the spectrum. For spatially con- 
tinuous systems supporting delocalized states the situ- 
ation is less obvious, but in case of a discrete energy 
spectrum a mapping onto a generalized tight-binding lat- 
tice applies as well for lower temperatures Ell ■ In 
any case, the configuration space variable can then be 
written as q(t) = a ■ s(t) with a typical length scale 
a and a dimensionless variable s(t) € {qi, . . . , q^\ with 
—S=q\ < <?2 < ■ ■ • < Id — S according to a d-level 
system (<fLS). Hence, the system Hamiltonian reads 



dLS 



= HE, — hS T 



(11) 



where E z describes the energetic distribution of the sites 
according to E z \q^) = e M | g M ) and S x the couplings be- 
tween them A M „ = (q^\S x \q^) , fi ^ v. In particular, in 
case of — q^ = 1 and nearest neighbor coupling only 
one recovers a (25 + l)-spin-boson model. 

For the discretization in time, the time axis is sliced 
via r uniformly spaced points with discretization steps 
t = t/r. The path integral in Eq. (JHJ) then becomes 



with 



P[{Sj}] 



J| K(s k ,s k+1 ) 



-HI*,}} 



(12) 



(13) 



The sum runs over all realizations of the discretized 
spin path {sj} = {s x = Sj, s 2 , . . . , s 2r , s 2r +i = Si}, and 
K(sj, Sj+x) denotes the coordinate representation of the 
free dLS propagation over the time interval r, i.e., 



K(s,s',t) = {s\exp(-iTH dhS /h)\s') 



(14) 



This propagator of the dLS Hamiltonian can be obtained 
from the eigenstates 



HdLs\4>a) = E a 



a 



as 



K(s,s', 



a=l 



S\(j>a)(<l>a\s' 



1, 



-irE a /h 



(15) 



(16) 



which can be easily computed numerically once the (iLS's 
parameters are specified. 

To arrive at a discretized form (in time) of the influence 
functional (J3J, the sum and difference coordinates 

V (t') = s(t')+s'(t'), Z( t ) = s(t')-s'(t') (17) 

are introduced, which read T)(t') = rjj = £j) for 

t' G [(J—1)t—t/2, (j — l)r+r/2] in their discretized form. 
The sum paths are also considered as "quasi-classical" , 
while the difference paths capture quantum fluctuations 
PJ. Equation © finally can be written as 



(18) 



i=2 



j>k=2 

where the kernels X^ Sl \ Xj-k, an d ^j-k follow from 
discretizing the twice integrated bath autocorrelation 
function Q(t) defined by Q(t) = L(t), Q(0) = with 
Q(0) = ifi/2. In the sequel we consider a spectral den- 
sity of the form 



(19) 



which is equivalent to ohmic damping with a cut-off fre- 
quency u> c . In this case Q(t) can be calculated analyti- 
cally and one obtains 



Q(t) - 2a 



ln(l + iu> c t) — In 



r(fi + it/hp)T(Q - it/tip) 



(20) 

with fi = 1 + l/(hf3uj c ) and the Gamma function T(z). 

Equations ^12(1. (|13|l and i|18|l constitute a discretized 
form of the populations JHl and thus provide a starting 
point for PIMC simulations. As it is well-known this 
method is handicapped by the dynamical sign problem 
|24|. It originates from quantum interferences between 
different system paths {sj}, causing a small signal-to- 
noise ratio of the stochastic averaging procedure. 

One approach to deal with this problem is based on the 
observation that the quasi-classical paths {rjj} in Eq. H18|) 
can be integrated out as a series of r — 1 matrix multipli- 
cations [2{| . This reduces the degrees of freedom from the 
2r— 1 variables {jj2<i<r+i> ^2<j<r} to the r — 1 quantum 
variables {^2<j<r} and therefore significantly improves 
the numerical stability of the corresponding MC simula- 
tions. While this technique, which in fact is yet another 
example of a blocking approach 26], works greatly for 
dissipative two- and three-state systems 20J, however, 
the increasing size of the corresponding matrices with 
the number of electronic sites turns these multiplications 
into a quite time consuming task. Since they have to be 
performed for every single update of the MC trajectory, 
the investigation of larger systems again requires excep- 
tionally long CPU times. 
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Nevertheless, this severe computational bottleneck can 
be profoundly alleviated on physical grounds 0, 12^ . 
Upon closer inspection one finds that the possibility of 
rewriting the integration over the quasi-classical coordi- 
nates in terms of simple matrix multiplications is due to 
the fact that the real-valued part of the bath autocorre- 
lation function JSJ governs only the quantum coordinates 
but not the quasi-classical ones. This real-valued part, 
which eventually leads to a damping out of quantum co- 
herences, introduces a non-local self-coupling among the 
quantum coordinates and is thus directly related to re- 
tardation effects, a main complication for treating dis- 
sipative quantum systems. Accordingly, the retardation 
effects influence the evolution of the quasi-classical coor- 
dinates contained in the imaginary part of the influence 
functional to a significantly weaker extend than that of 
the quantum coordinates. Neglecting them while gener- 
ating the MC trajectories leads to an only minor im- 
pairment of the sampling statistics, but causes an al- 
most complete decoupling between quantum and quasi- 
classical coordinates. This in turn allows for an enormous 
speed-up of the matrix multiplications. Accordingly, the 
sampling process could so far be accelerated by a factor of 
approximately 100 with respect to the ori gina l approach 
20] (for further details, we refer to Refs. |2l|, 22]), thus 
opening the door to treat the reduced dynamics of much 
larger and even many-body systems over sufficiently long 
times. 



IV. TIME-LOCAL MASTER EQUATIONS 

From the exact expression JHJ) for the reduced den- 
sity simplifications can be derived in certain limits in 
terms of time-local master equations. In the context dis- 
cussed here, these are of crucial importance, mainly since 
(i) master equations allow to study ranges in parameter 
space where the Monte Carlo sampling is rather expen- 
sive, e.g. for very weak or very strong friction, and (ii) 
they provide a basis to access the physical processes un- 
derlying the numerical data also by means of analytical 
techniques. 

In the weak friction limit one imposes that the level 
broadening due to friction is much smaller than k^,T and 
the typical level separation. It is thus natural to work in 
the basis spanned by the energy eigenstates of Hg and 
to treat the coupling Hj perturbatively. With increasing 
coupling, however, the environment drives the system to 
its pointer basis in which the system-bath coupling op- 
erator q is diagonal. As we have seen above in Sec. [HI 
in this basis the path integral formulation allows for a 
non-perturbative elimination of the reservoir degrees of 
freedom. In case that friction is very strong and thus 
level broadening much larger than level spacing and tem- 
perature, a master equation complementary to the weak 
friction range, the quantum Smoluchowski equation, is 
again available 0, 123 • 

What's about the intermediate regime? For TBSs dis- 



cussed here, substantial findings have been gained in the 
past decade 1] . Namely, it was shown that an exact re- 
tarded master equation exists, i.e., 

tdt't^t-t^p^it'), (2i) 

at „=l J 

where the kernels obey = — Ylv^u^i* 1 " Basically 
they represent a power series in the couplings A MI/ with 
corresponding spin-path integrals. To make the latter 
tractable, one applies the NIBA ^} or its generaliza- 
tion, the non-interacting cluster approximation (NICA) 
0, which has been shown to be accurate in a variety 
of parameter ranges. For instance, in case of an ohmic 
spectral density with a high frequency cut-off, it captures 
quantum coherence as well as pure relaxation dynamics. 
The corresponding kernels T ^ v read in lowest order in the 
intersite couplings 

f$(t) = 2A^ exp[-(^- 9;y ) 2 Q'(i)] 

cos [(e M - e v )t + - q v f Q"(t)] (22) 

with Q' = KeQ(t) and Q" = lmQ(t). The practical use 
of (|21|) is limited though due to the retardation. 

Now, for sufficiently strong dissipation and fast enough 
bath modes the kernel falls off on a time scale much 
shorter than the time scale on which the relevant reduced 
dynamics occurs so that we may set in (|21|l P V fli (t') ~ 
P v ,m(t) and 

IV = / dtt^it). (23) 
Jo 

Accordingly, (|2*T)l reduces to a simple Markovian rate 
equation 

P(i)=AP(t), (24) 

with a rate matrix A consisting of the individual rates 
T^ v . For nearest neighbor coupling, the golden rule 

rates T^Jl describe a sequential hopping process, while 
all higher order contributions to r^„ capture long-range 
hopping termed superexchange. Note that in contrast 
to weak coupling master equations as e.g. the Redfield 
equations, 1)2 1|) and l)24|) together with the corresponding 
transition rates apply also to the range of moderate to 
strong bath coupling and very low temperatures. 

As an explicit example of (|24|l we consider a one- 
dimensional tight-binding lattice with d = 7 sites, spac- 
ing 1, and constant nearest neighbor coupling corre- 
sponding to a spin-3-boson system (for details see [221]). 
This model can be seen as a simple realization of electron 
transfer alon g a molecular chain subject to a dissipative 
environment [2l], [28| . Initially, the charge is localized at 
the donor site Sj = q\ = — 3 and P s , _,g(t) monitors the 
dynamics towards the acceptor at qi = 3. Specifically, 
donor and acceptor are linked by a bridge with an impu- 
rity at its center. Donor/acceptor have vanishing onsite 
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FIG. 1: Populations along a molecular chain with d = 7 sites 
obtained from PIMC simulations (symbols) and from the local 
master equation approach (lines). Donor (D) and acceptor 
(A) are connected by a bridge with an impurity (I) at its 
center. The impurity has an onsite energy relative to D/A of 
ei/A = +5 (PIMC: diamonds; master equation: solid lines) 
and ei/A = —5 (PIMC: circles; master equation: diamonds). 
Other parameters are a = 0.1, Ahfi = 0.1, uj c / A — 5, see text 
for details. 



energy, the bridge sites are elevated by e#/A = 2.5, and 
the impurity site has ej/A = ±5. As seen in fig. 2] the 
exact PIMC data are very accurately described by the 
time-local master equation (|24[) despite the fact that all 
parameters are chosen such that one is close to a coher- 
ent /incoherent transition (for the given values of lj c and 
a coherences appear at inverse temperatures Tif3A 0.3 
and larger) and lo c /A — 5 is far from the scaling limit. 



V. DYNAMICS OF SPATIALLY CONTINUOUS 
SYSTEMS 

In order to go beyond the spin-boson type of models 
(one-dimensional tight binding lattice, equidistant spac- 
ing, constant nearest coupling), we start by presenting an 
approach to capture the dynamics of spatially continu- 
ous systems within the PIMC procedure outlined above. 
It turns out that this method applies to systems with a 
discrete energy spectrum. The basic idea is simple: For 
sufficiently low temperatures only the lowest lying states 
of the isolated system can be assumed to take part in the 
dynamics; thus, the full Hilbert space H$ can effectively 
be truncated to a subspace Hg spanned by the N low- 
est lying eigenstates. One then has to find a proper basis 
in this subspace, obtained by a unitary transformation 
from the eigenstate basis, in which a stochastic sampling 
of the path integrals (JSJ) can be performed. 

This sort of reduction is not new. In fact, the spin- 
boson model can be seen as originating from a double 
well potential, where only the two lowest lying states are 
taken into account. Of course, the goal here is to go be- 
yond: We want to capture the dissipative dynamics from 
very low up to moderate temperatures. For very high 
temperatures semi-classical or classical methods apply 
anyway. Further, we are interested in the regime of mod- 
erate friction, where neither of the known approximate 
formulations work so that a numerical treatment has to 
start from the exact reduced dynamics (JSJ. The proper 
basis for the sampling in the restricted Hilbert space is 
then the basis that diagonalizes the operator q in Hg K 
This idea to treat spatially continuous systems has been 
first applied in the numerical QUAPI approach 01 an d 
has been recently analyzed in [23j to derive a generalized 
master equation of the type given above 1|21|) . Hence, we 
sketch here only the central results briefly and refer to 
these previous works for further details. 



In the TV-dimensional subspace H s 
is defined as 



(AT) 



the pointer basis 



aqu \q v ) 



1, 



,N 



(25) 



and obtained from the eigenstate basis {\n}} by diago- 
nalizing a matrix with entries (n\q\m), n,m = 1, . . . , N. 
Its eigenvalues provide the dimensionless q v , while the 
eigenvectors are given as \q v ) — ~Y^b vn \n). The ba- 
sis set {| gv)} is also known as the DVR-basis (Dis- 
crete Value Representation). By representing the Hamil- 
tonian in this DVR-basis, one obtains "onsite ener- 
gies" he v = {q v \Hs\q v ) and "intersite-couplings" A„ M = 
{qv\Hs\qfj) /h. The original system is thus mapped onto 
a generalized A^-dimensional tight-binding lattice with 
non- equidistant sites at q v ,v = 1, . . . , N and non-nearest 
neighbor couplings A„ M . Now, to implement this TBS 
into the PIMC algorithm, one has to take into account 
the non-equidistant lattice by re-defining the £j and rjj 
variables properly. 

A rough estimate for the validity of the truncation pro- 
cedure to the lowest N eigenstates of the isolated system 
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FIG. 2: DVR-site populations for a harmonic oscillator with 
N = 5 eigenstates obtained from PIMC simulations (symbols) 
and a master equation approach (solid lines). Top left: Pi 
(diamonds) and P5 (circles), top right: P2 (diamonds) and P4 
(circles), bottom left: P3, bottom right: average position, see 
text for details. 



is given by the conditions that (i) the level broadening 
due to friction is of the same order as the level spacing 
or smaller and that (ii) the temperature is sufficiently 
low Nh(3uJo ^s> 1. While the second condition is obvious, 
the first one originates from the fact that for very strong 
friction the system tends to the classical limit again, see 

By way of example, we consider in the sequel a har- 
monic oscillator with mass M and freq uency ui p, so that 
the typical length scale is a = qo = y/h/MujQ. While a 
detailed account of the implementation into the PIMC 
algorithm and results of simulations for specific observ- 
ables, particularly in comparison with analytical results, 
will be given elsewhere, here, our main interest is to ex- 
plicitly show that apart from the conditions for the ap- 
plicability of time local master equations known for spin- 
boson models, the treatment of continuous systems via 
the described mapping comes with an additional com- 
plication. Formally, this issue has also been addressed 
recently in |23[ . Here, we directly compare numerically 
exact PIMC data with results from the simplified master 
equation. For this purpose the initial state is chosen as 
one of the DVR-states which suffices to study the time 
evolution of the populations on the sites q v . 

Results for N = 5 and N = 7 are depicted in figs. [3 El 
for a bath with oj c /ojq = 5, a = 0.2, and h0wp = 0.2. For 
this parameter set we have shown recently |2l], [2^ that 
time local master equations capture the exact dynamics 
of linear chains rather accurate even though one is close 
to an incoherent/coherent transition for the reduced dy- 
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FIG. 3: Same as in fig. © but for N = 7. Top left: Pi 
(diamonds) and P7 (circles), top right: P2 (diamonds) and Pq 
(circles), middle left: P3 (diamonds) and P5 (circles), middle 
left: P4, bottom: average position, see text for details. 



namics. The same is true here for N — 5, where we used 
(PUl together with the hopping rates it2*5|) . In contrast, 
for N — 7 deviations appear in the short to interme- 
diate time regime. For times of the order of 1/lu c and 
shorter (uot < 0.2), these are caused by adiabatic effects 
in the bath well-known from the dynamics in ordinary 
tight binding lattices. In the intermediate time range, 
however, where for certain populations pronounced max- 
ima occur, deviations must be ascribed a specific prop- 
erty related to the mapping onto a tight binding lattice 
with non-equidistant spacing. As mentioned above, this 
mapping requires a re-definition of the £ and rj variables, 
originally defined for a tight binding lattice with equidis- 
tant spacing of length a. 

Qualitatively, the situation is the following: Since the 
variance of the iVth eigenstate is roughly (Ag 2 )^ w Nq^ : 
the box in position space covered by N states has a typi- 
cal width L rs 2yNqo. The mean spacing between adja- 
cent lattice sites is thus a e g ~ L/N = 2qo/*/N. Hence, 
compared to a lattice with equidistant spacing of length 
a = qo, the re-definition of the £ and ry variables in the 
influence functional (|18fl comes with an additional factor 
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of the order (2/^/N) 2 , which basically renormalizes the 
damping kernel. Consequently, for fixed bath parameters 
and increasing TV, the effective coupling constant between 
bath and discrete system is not constant, but decreases as 
Aa/N . The system is thus driven into the range of weak 
coupling to the bath [23, where retardation effects due 
to quantum coherences become more substantial. Hence, 
the limit to a spatially continuous system N — > oo is non- 
trivial and must be performed by keeping a/N constant. 
Specifically, by comparing the spectral density 119(1 with 
the one usually introduced for spatially continuous sys- 
tems [l|, one finds 



the two-charge Hamilton operator reads 



4«7T 



7 
2uj 



(26) 



where 7 denotes the macroscopic damping constant ap- 
pearing in the classical Langevin equation. Note that the 
above relation cannot be seen as a strict equality since the 
spacing between adjacent DVR-sites varies slightly and is 
typically smaller deep inside the potential well. The value 
a = 0.2 chosen above corresponds for N = 5 to r y/u>o « 1 
and for N — 7 to 7/wo ~ 0.7. Now, a rough estimate for 
the validity of time local master equations can be gained 
by assuming that the exponential in (|22|l must fall off on 
a time scale sufficiently shorter than that of the dynamics 
of Pfj,(t). For moderate and low temperatures, this leads 
to l/(uj c ^a/N) < l/u>o and thus N ' < a(w c /cj ) 2 ■ In ac- 
cordance with this relation, the PIMC data presented in 
fig.[21ISlfor N = 5 can still be captured quantitatively by 
(|24jl . while for N = 7 only a qualitative agreement can 
be seen. To perform the continuum limit is thus not an 
easy task and certainly deserves further research. 



VI. CORRELATED TWO-PARTICLE 
DYNAMICS 



The dissipative real-time evolution of interacting many 
body systems has been left basically untouched so far. 
In certain limits, e.g. for very strong repulsive interac- 
tions, results could be derived from a simple hopping 
model p9| , but the intimate interplay between direct in- 
teraction, interaction mediated by the environment and 
dissipation has not been accessible. Here, for the first 
time we present PIMC results for the dynamics of two 
interacting particles along the lines described above. As 
in the previous section, in the sequel the formulation is 
outlined only briefly and we focus on conceptual prop- 
erties related to time-local master equations. Further, 
we consider the case of indistinguishable particles, called 
charges henceforth, so that the quantum nature of the 
particle statistics matters. Physically, the situation refers 
to interacting spinless fermions or interacting bosons. 

The free system is taken as a tight-binding lattice with 
spacing 1 and nearest neighbor coupling, where the two 
charges can be placed on d = 2S + 1 sites. Accordingly, 



H dLS - h 



where E z and Sx ,3 = a,b are straightforward general- 
izations of the corresponding operators introduced above 
for the single particle case, while U describes a site depen- 
dent Coulomb interaction specified below. The coupling 
to the bath is determined by the total dipol-moment of 
the system and thus follows directly from J2J: 



{E z - [£<■> + S^} + U) (27) 



H = H, 



(2) 



E 



p 2 

a 

2m r 



1 



2m a u> 2 



maljJt- 



(5*0 + Sf) 



(28) 



with Si f) \s^) = «Cfl j = a , b. 

Since the two charges are indistinguishable it is conve- 
nient for the PIMC simulation to work not in the single 
particle product basis, but rather in the basis of the many 
body states. For this purpose we introduce 



{| fl C«>)| s M>} — {| S) 3>} 



(29) 



with s < s. The advantage of this representation is that 

the d{d + l)/2 states {|s, §)} are orthogonal in contrast 

(2) 

to the original ones. Accordingly, ffl' g takes the form 



H 



(2) 



h [E z - S x + U] 



(30) 



where now 

E z \s,s) 
S x \s,s) 



U\s,s) 



(e s + e s )\s,s) 

S s=s (A s - 1 \s-l,s) + A s \s,s+l)) 
+ <5 s ^(A s _i|s - 1, s) + A s \s, s + 1) 
+ A s |s + l,s) +A 4 _i|M-l» 
-{u s ,s + ug, s )\s, s) . (31) 



For an onsite Coulomb energy this model is identical to 
a dissipative Hubbard model with two charges. Note, 
however, that for the PIMC simulations a non-local in- 
teraction can easily be taken into account. 

Now, as seen above, the free system dynamics enters in 
the discretized path integral formulation only through its 
short time propagator which is conveniently represented 
in the energy eigenbasis of ifVg in (|30|l . see (|16H . Fur- 
ther, since in (|28|l the two charges interact with the bath 
only via the sum S z a%> + S% the corresponding change 
to the many body basis in the influence functional is 
straightforward. It amounts to introduce discretized sum 
and difference paths according to 



(32) 
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FIG. 4: Mapping of a two charge system moving in a one- 
dimensional tight binding lattice onto a one particle system 
diffusing on a two-dimensional triangular lattice. 



such that the influence functional depends only on rj(t') + 
fj(t') and 

In principle, we could now start to implement the 
above representation into a PIMC algorithm. However, 
before doing so we go one step further and exploit the 
following crucial property: The dissipative dynamics of 
two indistinguishable charges on a one- dimensional lat- 
tice with d sites is equivalent to the dissipative dynam- 
ics of a single particle on a two-dimensional triangular 
lattice with d(d + l)/2 sites. To see this, one realizes 
that the coupling to the bath (si°^ + Si b ') /i£ with ±i£ 
as in Q can be written as ft ■ £ with the vector op- 
erators fi = (/i S* , fi Si^) and £ = (£,£), thus being 
identical to the system-bath coupling of a single parti- 
cle on a surface. Eventually, by a proper re-labeling of 
the populations P(si,Sf,Si,s~f;t) — > P(li,lf,t); U,lf = 
1, . . . , d(d + l)/2 one formally obtains the desired map- 
ping. This is illustrated in fig.0]where the actual coupling 
with the bath at a certain site is given by the projection 
of the corresponding "spin-vector" ft = (/js, /us) onto £ 
as just described and each site carries an onsite energy 
depending on e s + eg + (u s ,g + ug tS )/2. Due to the nearest 
neighbor-coupling, transitions can only occur in the ver- 
tical and the horizontal direction, respectively. Note that 
this mapping can be generalized to fermionic systems as 
well. 

To summarize, the advantages of the many body basis 
and the subsequent mapping onto a single particle 2d- 
lattice are: (i) the PIMC algorithm developed for single 
particles can be simply adapted to the case of two parti- 
cles, (ii) the configuration space to be sampled is reduced 
from d 2 to d(d + l)/2, and (iii) the master equations in- 
troduced in the previous section can be directly applied. 

Let us now analyze (iii) in more detail based on 
(|24[1 with the corresponding transition rates l|23|) . The 
simplest case is d = 2 with only three available sites 
(- 72, - 7 2 ), (- 72, 7 2 ), ( 72, 7 2 )- Obviously, for this 




4-4 c) 



FIG. 5: Dissipative dynamics of two correlated charges in 
a one-dimensional tight binding lattice with d = 5 sites ob- 
tained from PIMC simulations (symbols) and a master equa- 
tion approach (lines). Shown are the many-body state popu- 
lations P SyS i for the edge states of fig. @ with either s — —1.5 
or s' = 1.5. Circles, squares and diamonds denote situations 
where some sites have been removed from otherwise identical 
lattices as depicted in a), b), and c), respectively, referring 
to an increasing connectivity between the allowed sites. Solid 
(a), dashed (b), and dotted lines (c), respectively, depict the 
corresponding results from the master equation. 



case the triangular lattice coincides with a linear chain 
with three sites so that if a time local master equation 
applies for this latter case, it also applies for the two- 
charge case. For larger d, however, the situation changes 
fundamentally due to a different topology of the two- 
dimensional lattices which then contain bulk-sites with 
four adjacent sites (apart from d = 3 where there is no 
bulk site, but two edge-sites with three connections). Ac- 
cordingly, the typical dwell time on a bulk-site is con- 
siderably shorter than on a site with two connections, 
roughly, by a factor of 2. An environment which is able 
to destroy the phase coherence of the wave function on 
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each site of a linear chain before a jump to an adjacent 
site takes place, may be too slow to achieve the same on 
bulk-sites in a 2d-lattice. Hence, entangled many body 
states can survive on such a long time scale that a de- 
scription based on a master equation local in time fails. 
In fact, this can be seen in fig. 03 where PIMC data for 
d — 4 are depicted together with the corresponding dy- 
namics of the master equation l|24|l . By successively re- 
moving more and more bulk sites from the lattice, the 
agreement of the PIMC data with the prediction of the 
master equation, which is rather poor in case of the full 
lattice, increases, until upon creating a linear chain by 
removing all bulk states an almost perfect match is re- 
gained. 

The situation becomes even worse for a larger number 
of charges n since then the bulk sites of the n-dimensional 
cube attain 2n decay channels to adjacent sites, thus re- 
ducing the average dwell time by a factor of about 2n 
compared to the one-dimensional case n = 1. The con- 
clusion is the following: If for a bath characterized by 
w c and (3 a Markovian approximation applies for the dy- 
namics on a ld-lattice with nearest neighbor coupling 
A, i.e. fif3uj c < 1 and A/lo c <C 1, this approximation 
fails for the multi-charge dynamics, unless io c is taken 
to be very large and temperatures are sufficiently high, 
roughly, nA/uj c <C 1 and nh[3uj c < 1. Hence, for most 
cases only numerical approaches like the PIMC procedure 
presented here, seem to be applicable. 

VII. CONCLUSIONS 

The numerically exact PIMC approach has been 
pushed further to deal also with spatially continuous 
systems and correlated many-body dynamics. To bet- 



ter understand the numerical data in certain parameter 
ranges, to estimate the dissipative quantum dynamics 
before starting an involved PIMC calculation, and even 
to develop simplified models capturing the relevant pro- 
cesses, Markovian master equations are of great prac- 
tical use. To explore their applicability in the above 
situations, we had to consider generalized tight-binding 
lattices with either non-equidistant spacing or in higher 
dimensions. In both cases, additional restrictions must 
be imposed beyond the constraints known from equidis- 
tant one-dimensional TBSs so that for broader ranges of 
bath parameters no simple description seems to be pos- 
sible. However, it is surprising that in other domains 
comprising stronger dissipation and lower temperatures, 
Markovian master equations can still be found to work 
quite well, at least qualitatively. Thus, Markovian master 
equations together with numerically exact PIMC simula- 
tions, which often provide the sole mean to check for their 
validity in a certain scenario, provide a powerful means 
to study dynamical properties over the full time range. 
In particular, for the correlated many-body dynamics in 
presence of dissipation this may open the door to exam- 
ine the intimate interplay between Coulomb interaction, 
particle statistics, and phonon baths, e.g. for the charge 
transport in molecular and mesoscopic structures. 



Acknowledgments 

This work benefited from fruitful discussions with A. 
Komnik. We acknowledge financial support from the 
DFG through Grant No. AN336-1 and the Landesstiftung 
Baden- Wiirttemberg gGmbH. J. A. is a Hciscnbcrg fellow 
of the DFG. 



[1] U. Weiss, Quantum Dissipative Systems, Series in Mod- 
ern Condensed Matter Physics, Vol. 2 (World Scientific, 
Singapore, 1998). 

[2] C.W. Gardiner, Quantum Noise (Springer, Berlin, 1991). 

[3] G. Lindblad, Commun. Math. Phys. 48, 119 (1976). 

[4] A.G. Redfield, Adv. Magn. Reson. 1, 1 (1965). 

[5] J. Ankerhold, P. Pechukas, and H. Grabert, Phys. Rev. 
Lett. 87, 086802 (2001). 

[6] B. Vacchini, Phys. Rev. E 66, 027107 (2002); D. Baner- 
jee, B.C. Bag, S.K. Banik, D.S. Ray, Physica A 318, 6 
(2003). 

[7] J. Ankerhold, Europhys. Lett. 61, 301 (2003) ; J. Lucka, 
R. Rudnicki, and P. Hanggi, Physica A 351, 60 (2005). 

[8] J. Ankerhold, Europhys. Lett. 67, 280 (2004); L. 
Machura, M. Kostur, P. Hanggi, P. Talkner, and J. 
Lucka, Phys. Rev. E 70, 031107 (2004); J. Ankerhold, 
H. Grabert, and P. Pechukas, in: 100 years of Brownian 
motion, Chaos 15, 026106 (2005). 

[9] R.P. Feynman and F.L. Vernon, Ann. Phys. (N.Y.) 24, 
118 (1963). 

[10] A.O. Caldeira and A.J. Leggett, Physica 121A, 587 



(1983). 

[11] H. Grabert, P. Schramm, and G.-L. Ingold, Phys. Rep. 

168, 115 (1988). 
[12] R. Karrlein and H. Grabert, Phys. Rev. E 55, 153 (1997). 
[13] C.H. Mak and R. Egger, Adv. Chem. Phys. 93,39 (1996). 
[14] D. E. Makarov and N. Makri, Chem. Phys. Lett. 221, 

482 (1994). 

[15] W.T. Strunz, L. Diosi, and N. Gisin, Phys. Rev. Lett. 

82, 1801 (1999); J.T. Stockburger and H. Grabert, Phys. 

Rev. Lett. 88, 170407 (2002). 
[16] M.H. Beck, A. Jackie, G.A. Worth, and H.-D. Meyer, 

Phys. Rep. 324, 1 (2000). 
[17] A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P.A. 

Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 57, 

1 (1987). 

[18] A. Schmid, Phys. Rev. Lett. 51, 1506 (1983); M.P.A. 

Fisher and W. Zwerger, Phys. Rev. B 32, 6190 (1985). 
[19] R. Egger, C.H. Mak, and U. Weiss, Phys. Rev. E 50, 

R655 (1994). 

[20] R. Egger and C. H. Mak, J. Phys. Chem. 98, 9903 (1994); 
L. Muhlbacher and R. Egger, J. Chem. Phys. 118, 179 



10 



(2003). 

[21] L. Miihlbacher, J. Ankerhold, and C. Escher, J. Chem. 

Phys. 121, 12696 (2004). 
[22] L. Miihlbacher, J. Ankerhold, J. Chem. Phys. 122, 

184715 (2005). 

[23] M. Thorwart, M. Grifoni, and P. Hanggi, Phys. Rev. 
Lett. 85, 860 (2000); Ann. Phys. (New York) 293, 15 
(2001). 

[24] Quantum Monte Carlo Methods in Condensed Matter- 
Physics, M. Suzuki (ed.) (World Scientific, Singapore, 
1993), and references therein. 

[25] R. Egger and C. H. Mak, Phys. Rev. B 50, 15210 (1994). 



[26] C. H. Mak, R. Egger, and H. Weber-Gottschick, Phys. 

Rev. Lett. 81, 4533 (1998). 
[27] J. Ankerhold and H. Lehle, J. Chem. Phys. 120, 1436 

(2004). 

[28] J. Jortner and M. Bixon, eds., Adv. Chem. Phys. 106, 
107 (1999); A. Nitzan, Ann. Rev. Phys. Chem. 52, 681 
(2001). 

[29] E.G. Petrov and P. Hanggi, Phys. Rev. Lett. 86, 2862 
(2001); J. Lehmann, G.-L. Ingold, and P. Hanggi, Chem. 
Phys. 281, 199 (2002). 



